{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Curve Fitting\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "Try this Jupyter Notebook Modeling Example to learn how you can fit a function to a set of observations. We will formulate this regression problem as a linear programming problem using the Gurobi Python API and then solve it with the Gurobi Optimizer.\n",
    "\n",
    "This model is example 11 from the fifth edition of Model Building in Mathematical Programming, by H. Paul Williams on pages 266 and 319-320.\n",
    "\n",
    "This modeling example is at the beginner level, where we assume that you know Python and that you have some knowledge about building mathematical optimization models. The reader should also consult the \n",
    "[documentation](https://www.gurobi.com/resources/?category-filter=documentation) \n",
    "of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br />\n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$i \\in \\text{Observations}=\\{1, .. ,n\\}$.\n",
    "\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$x_{i} \\in \\mathbb{R}$: Independent variable value at observation $i$.\n",
    "\n",
    "$y_{i} \\in \\mathbb{R}$: Dependent variable value at observation $i$.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$a \\in \\mathbb{R}$: Value of the constant term in the function that explains the values of $y$ in terms of the values of $x$.\n",
    "\n",
    "$b \\in \\mathbb{R}$: Coefficient of the linear term in the function that explains the values of $y$ in terms of the values of $x$.\n",
    "\n",
    "$u_{i} \\in \\mathbb{R}^+$: Positive deviation of the proposed function of x with respect to the value of y at observation $i$.\n",
    "\n",
    "$v_{i} \\in \\mathbb{R}^+$: Negative deviation of the proposed function of x with respect to the value of y at observation $i$.\n",
    "\n",
    "$z$: Value of the maximum deviation.\n",
    "\n",
    "We model the problem for the first goal:\n",
    "\n",
    "* Fit a line $y=a+bx$ to the given data set in order to minimize the sum of absolute deviations of each observed value of $y$ from the value predicted by the linear relationship.\n",
    "\n",
    "\n",
    "### Constraints Problem 1\n",
    "\n",
    "**Deviation**: Each pair of corresponding data values $(x_{i},y_{i})$ gives rise to the following constraint.\n",
    "\n",
    "\\begin{equation}\n",
    "bx_{i} + a + u_{i} - v_{i} = y_{i}  \\quad \\forall i \\in \\text{Observations}\n",
    "\\end{equation}\n",
    "\n",
    "Where $x_{i}$  and $y_{i}$  are  the given values in the set of observations, $b$, $a$, $u_{i}$ and $v_{i}$ are variables. \n",
    "The positive deviation $u_{i}$ and the negative deviation $v_{i}$ give the amounts by which the values of $y_{i}$ proposed by the linear expression differ from the observed values.\n",
    "\n",
    "### Objective Function Problem 1\n",
    "\n",
    "**Total deviation**: The objective is to minimize the total positive and negative deviations.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Minimize} \\quad \\sum_{i \\in \\text{Observations}} (u_{i} + v_{i})\n",
    "\\end{equation}\n",
    "\n",
    "We now provide a model formulation for the second goal:\n",
    "\n",
    "* Fit a line $y=a+bx$ to the given data set in order to  minimize the maximum deviation of all the observed values of $y$ from the value predicted by the linear relationship.\n",
    "\n",
    "For this new formulation, in addition to the \"Deviation constraints\", we need to include the following constraints.\n",
    "\n",
    "### Constraints Problem 2\n",
    "\n",
    "**Maximum deviation**: The following constraints ensure that the decision variable $z$ takes the value of the maximum deviation.\n",
    "\n",
    "\\begin{equation}\n",
    "z \\geq u_{i}  \\quad \\forall i \\in \\text{Observations}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "z \\geq v_{i}  \\quad \\forall i \\in \\text{Observations}\n",
    "\\end{equation}\n",
    "\n",
    "### Objective Function Problem 2\n",
    "\n",
    "**Minimum/Maximum deviation**: The objective is to minimize the maximum deviation.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Minimize} \\quad z\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.7.0 & Gurobi 9.1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data\n",
    "\n",
    "We define the corresponding values for $x$ and $y$ in the set of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample data: values of independent variable x and dependent variable y\n",
    "\n",
    "observations, x, y = gp.multidict({\n",
    "    ('1'): [0,1],\n",
    "    ('2'): [0.5,0.9],\n",
    "    ('3'): [1,0.7],\n",
    "    ('4'): [1.5,1.5],\n",
    "    ('5'): [1.9,2],\n",
    "    ('6'): [2.5,2.4],\n",
    "    ('7'): [3,3.2],\n",
    "    ('8'): [3.5,2],\n",
    "    ('9'): [4,2.7],\n",
    "    ('10'): [4.5,3.5],\n",
    "    ('11'): [5,1],\n",
    "    ('12'): [5.5,4],\n",
    "    ('13'): [6,3.6],\n",
    "    ('14'): [6.6,2.7],\n",
    "    ('15'): [7,5.7],\n",
    "    ('16'): [7.6,4.6],\n",
    "    ('17'): [8.5,6],\n",
    "    ('18'): [9,6.8],\n",
    "    ('19'): [10,7.3]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. The variables of the model are the constant term and coefficient of the linear term of the function f(x), the positive and negative deviations, and the maximum deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('CurveFitting')\n",
    "\n",
    "# Constant term of the function f(x). This is a free continuous variable that can take positive and negative values. \n",
    "a = model.addVar(lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=\"a\")\n",
    "\n",
    "# Coefficient of the linear term of the function f(x). This is a free continuous variable that can take positive \n",
    "# and negative values.\n",
    "b = model.addVar(lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=\"b\")\n",
    "\n",
    "# Non-negative continuous variables that capture the positive deviations\n",
    "u = model.addVars(observations, vtype=GRB.CONTINUOUS, name=\"u\")\n",
    "\n",
    "# Non-negative continuous variables that capture the negative deviations\n",
    "v = model.addVars(observations, vtype=GRB.CONTINUOUS, name=\"v\")\n",
    "\n",
    "# Non-negative continuous variables that capture the value of the maximum deviation\n",
    "z = model.addVar(vtype=GRB.CONTINUOUS, name=\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each pair of corresponding data values $x_{i}$ and $y_{i}$ gives rise to a constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deviation constraints\n",
    "\n",
    "deviations = model.addConstrs( (b*x[i] + a + u[i] - v[i] == y[i] for i in observations), name='deviations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective function of problem 1 is to minimize the total positive and negative deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function of problem 1\n",
    "\n",
    "model.setObjective(u.sum('*') + v.sum('*'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 19 rows, 41 columns and 75 nonzeros\n",
      "Model fingerprint: 0x0bec2f7b\n",
      "Coefficient statistics:\n",
      "  Matrix range     [5e-01, 1e+01]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [7e-01, 7e+00]\n",
      "Presolve removed 0 rows and 1 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 19 rows, 40 columns, 75 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0      handle free variables                          0s\n",
      "      20    1.1466250e+01   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 20 iterations and 0.01 seconds\n",
      "Optimal objective  1.146625000e+01\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('CurveFitting.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "_________________________________________________________________________________\n",
      "The best straight line that minimizes the absolute value of the deviations is:\n",
      "_________________________________________________________________________________\n",
      "y = 0.6375x + (0.5813)\n"
     ]
    }
   ],
   "source": [
    "# Output report\n",
    "\n",
    "print(\"\\n\\n_________________________________________________________________________________\")\n",
    "print(f\"The best straight line that minimizes the absolute value of the deviations is:\")\n",
    "print(\"_________________________________________________________________________________\")\n",
    "print(f\"y = {b.x:.4f}x + ({a.x:.4f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Problem 2, it is necessary to introduce another variable $z$ to capture the value of the maximum deviations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Maximum deviation constraints\n",
    "\n",
    "maxPositive_deviation = model.addConstrs( (z >= u[i] for i in observations), name='maxPositive_deviation')\n",
    "\n",
    "maxNegative_deviation = model.addConstrs( (z >= v[i] for i in observations), name='maxNegative_deviation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective function for Problem 2 is to minimize the maximum deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 57 rows, 41 columns and 151 nonzeros\n",
      "Coefficient statistics:\n",
      "  Matrix range     [5e-01, 1e+01]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [7e-01, 7e+00]\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    0.0000000e+00   1.146625e+01   0.000000e+00      0s\n",
      "      11    1.7250000e+00   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 11 iterations and 0.01 seconds\n",
      "Optimal objective  1.725000000e+00\n"
     ]
    }
   ],
   "source": [
    "# Objective function for Problem 2\n",
    "\n",
    "model.setObjective(z)\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "_________________________________________________________________________________\n",
      "The best straight line that minimizes the maximum deviation is:\n",
      "_________________________________________________________________________________\n",
      "y = 0.6250x + (-0.4000)\n"
     ]
    }
   ],
   "source": [
    "# Output report\n",
    "\n",
    "print(\"\\n\\n_________________________________________________________________________________\")\n",
    "print(f\"The best straight line that minimizes the maximum deviation is:\")\n",
    "print(\"_________________________________________________________________________________\")\n",
    "print(f\"y = {b.x:.4f}x + ({a.x:.4f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
